//calc_row_pvalue.c
//Programmer: Rory Carmichael
//Date: 2/17/2011
//Purpose: row_wise calculation executable for masterworker

#include <stdio.h>
#include <unistd.h>
#include <ctype.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <math.h>

float* get_peak_frequencies(char *filename, int *numcols) {
	FILE *infile;
	int colsize = 0;
	float *pvals= NULL;
	int i = 0;
	colsize = *numcols;
	infile = fopen(filename, "rb");
	if ( infile != NULL ) {
		fread(&colsize,sizeof(int),1,infile);
		pvals = (float *) malloc(colsize*sizeof(float));
		fread(pvals,sizeof(float),colsize,infile);
		*numcols = colsize;
		return pvals;
	}
	return NULL;
}

int is_set(int position, uint32_t n) {
	return n & 1<<position;
}

int compare_float(const void *a, const void *b) {
	float fa = *(float *) a;
	float fb = *(float *) b;
	return (int) (fa - fb);
}

int compare_int(const void *a, const void *b) {
	int ia = *(int *) a;
	int ib = *(int *) b;
	return (int) (ia - ib);
}

float *calc_cell_pvalue(int dist, float *pvallist, int numcols) {
	int i = 0;
	float *pcalc = NULL;
	pcalc = (float *) malloc(numcols*sizeof(float));
	for(i=0;i<numcols;i++) {
		pcalc[i] = 1 - pow( (1-pvallist[i]) , (2*dist)+1 );
	}
	return pcalc;
}

char *as_binary(uint32_t num) {
	char *bitstring;
	bitstring = (char *) malloc(33*sizeof(char));
	bitstring[32] = '\0';
	int i = 0;
	while(i < 32) {
		if (num & 1<<i) {
			bitstring[31-i] = '1';
		} else {
			bitstring[31-i] = '0';
		}
		i++;
	}
	return bitstring;
}

float bernoulli_exact(int tgt_num, float *pcalc, int pcalc_len) {
	int sum = 0;
	float microprob = 1.0;
	float semiprob = 0.0;	
	float bitprob;
	uint32_t max = (uint32_t) pow(2, pcalc_len);
	uint32_t i;
	int j;
	for(i=0; i<max; i++) {
		sum = 0;
		microprob = 1.0;
		sum = __builtin_popcount(i);
		for(j=0; j < pcalc_len; j++) {
			if(is_set(j, i)) {
				bitprob = pcalc[j];
			} else {
				bitprob = 1 - pcalc[j];
			}
			microprob *= bitprob;
		}
		if(sum >= tgt_num) {
			semiprob += microprob;
		}
	}
	if(semiprob == 0) { semiprob = 1; }
	return semiprob;
}

float calc_row_pvalue(int numcols, int *distlist, float *pvallist) {
	int i = 0;
	int j = 0;
	float masterprob = 1.0;
	float semiprob = 1.0;
	float *cellprobs;
	qsort(distlist,numcols,sizeof(int),compare_int); 
	for(i = 0; i < numcols; i++) {
		cellprobs = calc_cell_pvalue(distlist[i],pvallist,numcols);
		semiprob = bernoulli_exact(i,cellprobs,numcols);
		if(semiprob < masterprob) {
			masterprob = semiprob;
		}
	}
	return masterprob;
}

int main (int argc, char **argv) {
	FILE *infile;
	int i = 0;
	int j = 0;
	float *pvallist = NULL;
	char *filename = NULL;
	char *pvalfile = NULL;
	char *result = NULL;
	char delims[] = "\t\n";
	char line [ 10000 ];
	int *distlist = NULL;
	int *curdists = NULL;
	float *curpvals = NULL;
	int numcols = 0;
	int gapcount = 0;
	float pval = 1.0;
	if (argc != 3) {
		fprintf(stderr, "USAGE: calc_row_pvalue pvaluefile linefile");
		exit(1);
	}
	pvalfile = argv[1];
	filename = argv[2];
	pvallist = get_peak_frequencies(pvalfile, &numcols);

//	printf("PVALLIST: %f",pvallist[0]);
//	for(i=1;i<numcols;i++) {
//		printf("\t%f",pvallist[i]);
//	}
//	printf("\n");

	distlist = (int *) malloc(numcols*sizeof(int));
	if(pvallist != NULL) {
		infile = fopen(filename, "rt");
		if(infile == NULL) {
			fprintf(stderr, "Failed to open %s\n", filename);
			exit(1);
		}
		while (fgets(line,sizeof(line),infile) != NULL ) {
			char ugly [ 10001 ];
			sprintf(ugly,"%sx",line);
			if(strcmp(ugly,"x") == 0) {
				break;
			}
			result = strtok(line,delims);
			i = 0;
			gapcount = 0;
			while (result != NULL && (strcmp(result, "\n") != 0) && (strcmp(result, "") != 0)) {
				if (i > numcols) {
					fprintf(stderr, "MALFORMED INPUT FILE... %d (with value %s)is too many columns\n", i, result);
					exit(1);
				}
				if (strcmp(result, "-") != 0) {
					sscanf(result, "%d", &distlist[i]);
				} else {
					distlist[i] = -1;
					gapcount++;
				}
				i++;
				result = strtok(NULL, delims);
			}
			curdists = (int *) malloc((numcols-gapcount)*sizeof(int));
			curpvals = (float *) malloc((numcols-gapcount)*sizeof(float));
			i = 0;
			j = 0;
			while(i < numcols-gapcount) {
				if (distlist[j] != -1) {
					curdists[i] = distlist[j];
					curpvals[i] = pvallist[j];
					i++;
				} 
				j++;
			}
			i = 0;
			pval = calc_row_pvalue(numcols-gapcount, curdists, curpvals);
			fprintf(stdout, "%f\n", pval);
			free(curdists);
			free(curpvals);
			curdists = NULL;
			curpvals = NULL;
		}
	}
	return 0;
}
